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Abstract 


Closed-loop, analytical solution of Linear, algebraic equation containing many unknown variables are found in many 
mathematical modeling equations and network analysis involved in healthcare and neuroscience research. Finding 
unique, analytical solution for linear, algebraic equations has diverse application in many fields. Computation of 
unique, non-trivial solutions to large array containing several variables is a computationally overwhelming task. The 
paper begins by introducing the concept of network analysis in modelling for healthcare and neuroscience research. 


A simple example of network modelling is the logistics of delivering vaccine like COVID19 from a company to . 


Primary health center while maintaining cold-storage of the vaccine. This is followed by explanation of the working 
of computationally efficient, best-practice LU Decomposition with partial pivoting algorithm to solve dense linear 
equations. This is followed by narration of building, testing, obfuscation, compiling and release of an Android Graphic 
user interface application implementing the above methods. 

The final part of the paper examines the exceptional accuracy and efficiency of solving, dense matrix equations on 
Android Run Time machine using this approach. The calculated Poisson modelling of probability of Stochastic, 
Singularity event is = 3.26 x 10~°per floating point operation a number derived after running 15,625, 000,000 
floating-point operation runs. The mean execution time was 88.534 seconds , for solving matrix equation 
[A]withN=60 variables in performing N?=216000 computations. The whole working Android application containing 
many other tools is hosted on the GitHub and Figshare platforms along with additional graphs, dataset, Java, and 
Python programs used to complete this study. 


KeyWords:LU Decomposition, Matrix algebra, Mathematical modelling, Computer aided algebra, Linear algebra, 
Linear equations 
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Introduction There are numerous situations in science which 


First degree equations of two to any number of require analytical, solutions a system of Linear 
unknown variables, all of which have an exponent equations.There are numerous instances in diverse 
of 1 (x1 so the highest degree is 1) are called Linear fields like Ecology, Plant taxonomy, Business 
systems of equations. The word “Linear” here may mathematics, Biostatistics, Precalculus, Electronics: 
remind the reader that these equations are Kirchhoff 's current law, Physics: 

represented by straight lines (or planes) in 

orthogonal coordinate systems. 
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Simple Mass-Spring Systems, Linear Programming, 
linear econometric models, mathematical 
modelling, solutions of elliptic partial differential 
equations and so on. In addition, even many non- 
linear systems can be modelled with a reductionist 
approach to system of Linear equations, for 
example Numerical or approximate solutions of 
partial Differential equations. They also appear in 
digital image processing and compression, digital 
signal processing, digital noise filtering in systems 
and operational theory of analysis as stated by 
Yang, M.-H., Yang, F.-Y. and Oyang, Y.-J. (2013). 
Linear equations and Network analysis in 
Healthcare and neuroscience research 
Finding unique solutions to Linear, algebraic 
equations are a recurring theme in many 
mathematical modelling and network analysis 
problems in healthcare and neuroscience research. 
For example, logistics supply chain model of 
COVID19 vaccination is modelled in figure1. The 
supply chain and warehousing of the many vaccines 
including Coronavirus is difficult problem, as needs 
continuous maintenance of cold _ storage 
temperature in the warehouse called the Cold 
chain. Coronavirus vaccines like the Indian Covaxin 
and Covishield need to be stored 2-8 degrees, and 
by contrast Pfizer’s vaccine needs to be stored at - 
70°C. Managing logistics of “cold chain” is crucial to 
success of vaccination campaign. This simple 
example can be network modelled using a rooted, 
directed tree graph as shown below in figure1. This 
graph was created using free Graphing software 
Gephi as stated by Bastian, M., Heymann, S. 
&Jacomy, M., (2009). This is also a directed acyclic 
graph (DAC) where the edges are directed but there 
is no rotation. 
GraphGyaccineLogistics = (V,E)N, = IV] =7, Ne 

= |E|=8 
In the Figure1. simple example biological network 
Graph of order (nodes) 7 and Size (edges) 8. This 
network is defined by nodes or vertices are number 
vaccine available at each center utilized or wasted. 
The vaccine availability and flow are independent 
of the nodes. Since the chance of breaking the cold 
chain increases with increasing number of nodes 
(by repeated storage and mobilization in and out of 
freezer of for dispensing) the shortest walk or path 
through this example tree called the Geodesic 
would be the optimum logistical solution. In 
contrast, the longest path called Diameter would be 


the most logistical critical path in this graph. Such a 
healthcare network can be _ mathematically 
modelled with the following sparse matrix 
equation, called the Adjacency matrix which is 
basically the summary of fundamental connectivity 
within the network as described by Kolaczyk, E.D. 
(2010). 
Adjacency matrix of 

1ifi,j € (belongsto) Edgeset 


[A]a,; = 


Oifi,j € (doesnotbelongthe) Edgeset 
Another useful matrix is Incidence matrix which 
contains directional information about the edges of 
the network. Incidence matrix is defined for a 
directed graph as following, based on paper of 
Seshu S. and Balabanian, N. (1967). 


Incidence matrix[A]a;; = 


1 ifi,j € (belong)Edgeand1 ifarrowawayfromvertexa 


=I 

0 ifi,j € (notbelong)Edgeset 
Thus, the Incidence matrix [/n] is a square matrix, 
computed for graph [A] like this based on the above 
rule. 


Node1;1 -1 0O 0 0 0 0 0 
Node2 |-1 1 1 1 0 0 0 0 
Node3a| 0 0 0 -1 O O +1 +1 
Node4a| 0 0 0 0 0 0 -1 O 
Node3b| 0 0 -1 0 41 O 0 -1 
Node4b| 0 0 0 0 -1 +1 O 0 
Node5Sb' 0 0 0 0 0 -1 O 0 
Node1i;s1 -1 0 0 0 0 0 0 
Node2 |-1 1 1 1 0 0 0 0 
Node3a|0 0O O -1 0 0O +41 +1 
Node4a| 0 0 0 0 0 0 -1 O 
Node3b| 0 0 -1 0 41 0O 0 -1 
Node4Ab| 0 0 0 0 -1 +1 O 0 
Node5Sb' 0 0 0 0 0 -1 O 0 

c1 
c2 
c3 
=]c4 
c5 
c6 
c7 


ormatrix, condensedform|[A] x [X] = 


or —1ifarrowdirectedtowardsvertexa 
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Fig1. A directed, rooted tree graph showing simple network analysis of logistics in distribution of 
vaccines like Covid-19 from a district hospital node(2), sourced from production, the root node 
(edge 1) and returned to root node (edge 2) if unused, node(3a) is cold-storage district hospital 
(district walk-in freezers), and node (4a) are people vaccinated directly in the hospital and node 
are those vaccinated in the community (5b) after storage node (4b) Small freezers in the local 
Primary Health centers node (3b). Node. Primary Health centers occasionally gets stocks from 
district walk in freezers (edge 8). This healthcare situation can be mathematically modelled using 
Network analysis and Matrix equations, this is based on Khanna, V.V., Chadaga. K. et. al. (2022) 
and Allen, L.J.S. (2007) 


The first row is given by +X1(¢) — X2(¢) + 0X3) + 
Ox4(t) + OXs(t) + OX6(t) = C1 so on leading to N=8 
linear equations which must be analytically solved 
by Gaussian elimination or any other methods, and 
unique values of x4(¢),%2(¢), calculated in a method, 
as described by Ma'‘ayan A and Macarthur, B.D 
(2012).. This situation being aided by sparse nature 
of the matrix and the nodal relations are not so 
simple in many real-life situations. With the 
increasing number of variables, the number of 
steps required to manually solve this system of 
equation by for example Gauss Elimination (by 
reducing the matrix to reduced row echelon form 
by series of elementary row transformations, as 
detailed by Yang, W., Zhang, D. et.al. (2021)) 
becomes very large. Most of the manually solved 
examples in textbooks contain real numbers from 1 
to 10, which are hardly ever found in real life 
equations. 


Objectives 


Analytically solve dense linear equations of size (N) 
limited only by Operating system limitation or RAM 
on Android cell phones using LU Decomposition 
and Partial pivoting in Android Java language 
(Android studio). 


Methodology 


Generally, the unknown variables in the System of 
Linear systems of equations are represented by x 
(same as x, y, z or 3 variables in standard books). 
We designate a square matrix [A] containing N 
unknows of size N? where N = {n € Nandn> 2} 
Variables/Unknowns 
X41 (x), X2(y), X3(Z), X4, Xs «...Xywhere n can be any 
natural number (but not zero) 
¥ = (Gq), Ge), (og), (aa) ~ Gen) 
Column matrix [X] is also called a vector of 
dimension or axis one and variables on this axis are 
called the components of this vector. The number 
of N is called the dimension of the vector space as 
stated by Abadir, K.M. and Magnus, J.R. in their 
paper (2005). 
Common equations are thus written in this 
convenient form for the purpose of this paper and 
computer program. 

ax+by+dz=k which is synonymous 
with 44X41 + Q42X42 + Q143X3 =c,so on for ith 
equation thus ail + ai2 +... or Summation of 
pres = c; . The term Cn is a column matrix (1 
column only) also called a vector C thus Ic] and 
inar thus 


column matrix the equation is 
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represented by vector dot product iba [A] = Ici. 
This is called the Dot product Vector representation 
of a linear equation (as described at this site 
http://www.wolframalpha.com (1996). 

Matrix Representation of an equation is the 
depiction of a linear equation as a product of a 
square matrix with 2 column matrices as stated by 
Stewart. G.W. in his book (2000). 


Q41 42. A43] [%1 Cy 
[A] =|]@21 422 | [a = a 
Q23 432 4331 1X3 C3 
So, on... 
[aia Qi2 diz....a;;| [x;] = [cj] 
Or simply 
[A] - [X] = [C] 


Simply rearranged 
[X] = inverse([A]) x [X] = [A]71 x [C] 
adj[A] 

= Getlal * 1) 
On the precondition that [A] is a square non- 
singular matrix. 
For the sake of simplicity, the above notation is 
going to be used throughout this paper, with 
arrows omitted and is a synonymous with any 
other representation as stated by Xian-Da Zhang 
(2017). 
Another representation is Augmented matrix 
representation, this is not used as this is notation 
for hand solution by Gaussian elimination by 
reducing to row echelon form by performing hand, 
elementary row transformation on both sides of the 


matrix. 
X44 X4Ayq -X4A43| Cy 
[A] = |%2@21 %2A22 X2Q43 | | 
X3A23 X3A32 %3433| C3 


Diagonal Matrix is one in which all elements below 
the diagonal line are zero, but diagonal line itself is 
nonzero. An Identity matrix has 1 in the diagonal 
line. 


aj, O 0 1 0 0 
diag[A]=|0 a2 O |IdentityJ=|0 1 O 
0 0 a33 0 0 1 


A permutation matrix P is a binary matrix (of 1, 0 
only) which has 1 in each row and each column. 
This matrix (as we will see further) is an Operator 
with causes interchange of rows, like this example. 
? 1 I" b) _ B d 

1 0 c ad a b 
A sparse matrix is a type of matrix with sufficient 
zeros so that they can be manipulated in a 
computationally efficient manner. In contrast a 
dense matrix is neither is repetitive nor sparse to 
help in reaching a unique solution. The objective is 
to solve this dense system of linear equation in 
finite number of steps in a most efferent way called 
closed-loop solutions as described by Debonis, M.]. 
(2022). Open loop or trivial solution are of not 
much interest like for example +; = 0,4, =0,xy = 
0, will solve any equation (even quadratic 
equation) but are not of any real utility as stated by 
Bronson R. and Costa, B.G. (2009). 


Effect of a Permuation Matrix [P] 


aon ee 


aaa 


Rows ofa Matrix are Switched 
Fig.2. Schema of Partial Pivoting by the row-shifting effect of product with Permutation matrix [P] 


Methods of arriving at closed-loop solution of 


Matrix equation [x] x Al = Ic] 

There are many ways such a dense, matrix equation 
maybe solved which are broadly classified into: 

1. Stanimirovi, I. (2017) observed that Analytical 
methods or Direct methods by following a specified 
algorithm producing accurate solutions to [X] x 
[A] =[C] , these are useful for dense, non- 


sequential matrix equations, examples are Cramer’s 
rule, Inversion of Matrix, Gauss Jordan method of 
direct elimination by inspection, and LU 
Decomposition. 

2. Young, D.M. (2003) stated that Numerical 
methods or Iterative methods by sequence of 
successively better and better approximations of 
the solutions, these are useful for sparse, 
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sequential, or structural matrix like Toeplitz’s, 
Hankel’s, Vandermonde, Cauchy, Pick-ups, Bzout’s 
and Toeplitz’ matrices. To reminisce a few extant 
numerical methods: Gauss-Seidel method, SOR, 
SSOR, Richardson’s method, and Jacobi method. 


Problem with common classroom methods of 
solving large Linear equations 


Golub, G.H., Chan, R.H., Greif, C., O’Leary, D.P. and 
Ebrary, I. in their paper (2007) observed that for 
example, solution of Matrix equation of order n 
using rearranged [X] = inverse([A])[C], 
computation of [A]~+ of a square matrix of size N 
will need to solve (N+1) number of determinants 
and performN + 1! factorial number of operations 
in terms of multiplications and additions. For 
example, solving a linear equation having 6 terms 
will need 7! number of operations which will 5040 
accurate steps by hand. Analytical, closed-loop 
solution of linear equation by Cramer’s rule with 10 
terms will need a staggering 3,628,800 steps. 


Lower and matrix 


decomposition 


Upper _ Triangular 
LU decomposition is attributed Alan Turing in 
1948, who showed the factorization of a Matrix into 
2 triangular, Matrices: The Upper [U] and Lower 
[L]matrix is unique (proof is given in all textbooks) 


Upper 


LU Decomposition 
iilustrated 


for any square, non-singular matrix. LU 
factorization is the computationally most efficient 
algorithm for solving such a system of linear 
equations, in contrast to other methods in 
literature like Block LU decomposition, Rank 
factorization, Cholesky decomposition, QR 
decomposition, and singular value decomposition. 
For instance, Cholesky decomposition produces 
contain square roots which are computationally 
difficult to solve and QR decomposition method 


4], N®? takes twice the number of computations to 


achieve the same _ results as 
AnatolijDmitrievicMyskis (1972). 
LU Decomposition/ factorization is a more efficient 
implementation of Gauss direct elimination 
method. This efficiency allows for solving Matrix 
equations on a mobile processor, a task which was 
relegated to Supercomputers in the past. This is 
factorization scheme is shown in picture that 
follows. 

Any Square Matrix as factors of 2 Square matrices. 
Such an operation is possible provided matrix is 
non-singular. 

[A] = [L] x [U] 


observed by 


Such a factorization is particularly beneficial in 5 


solving revolving differential equations whose [C] 
direction of Vector changes due to revolutions of 
the vector. 


Triangular matrix 


Lower Triangular Matrix 


Fig.3. Schema of Factorizing a matrix [A] into an Upper and Lower triangular matrix for purpose of 
fast and efficient computation 


Special properties of diagonal elements of L and U 
matrices making computations easy 

Computing determinants is easy for Upper and 
lower Triangular matrices and involves Product of 
the Diagonal elements of the matrices. 


a4, 2 4 
|A| = 0 a22 7 |= Qay1 x a22 x a33 x Agg...X 
0 0 a3 


aij 
Similarly inverse of Lower Triangular matrix is 
simply the sign inversed elements below the 
diagonal. 

[L]-* = (-)) x [L] 
Solution of [X] is easily reduced to finding inverse 
of [U]~1 and matrix multiplications. 

[X] = [A]-* x [C] = [L]-* x [U}"* x [Cc] 
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LU Decomposition elegantly, simplifies the closed- 
loop solutions of complex, system of equations. So, 
the Upper and lower Triangular matrices can be 
easily used to find several other closed-loop 
solutions like Determinant of|A|, Inverse of matrix 
[A]~'so on. This flexibility, ease and computational 
efficiency makes the LU factorization the preferred 
analytical method of solving matrix equations as 
described by A. G. Kurosch (1988). 

Smaller number of floating-point operations (fops) 
or steps making computations easy 

The number of floating-point operations (fops) or 
steps need to complete LU factorization of a dense, 
Linear equation of order N is given by the formula 
O(N?).The number of floating-point operations 
(fops) needed to solve a matrix of size N is given by 
the equation O(N?)when W € can be any Natural 
number WN > 2. and NW —1 division operations, 
which are computationally complex. Therefore, an 
equation with 6 unknowns will need only 216 steps 
and equation with 10 unknowns can be solved in 
100 accurate steps, for 20 it is 4000 steps only as 
stated in paper by Bisseling, R.H. (2020). In 
contrast, the LU factorization process of large 
numerical arrays is not without problems as we 
shall see in next section. 

Computational problems inherent to factorization 
of a Matrix into U and L Matrices 

Non-Existence of a Lower and Upper Triangular 
matrix 

The LU-decomposition is computational efficient, 
best-case algorithm but this method may fail in 
some, peculiar matrices. 

Problem of Singulate or degenerate Matrix and 
Singularities in calculation 

A matrix [A] is non-singular on the condition that 
there exists such a matrix 


[A] x [A]~* = [4]-* x [4] = [1] 

[I] is the identity matrix, element 1 along the 
diagonal and rest all the non-diagonal elements are 
zero. A singular matrix is also called degenerate 
matrix and determinant is zero |A| =0, such a 
system of matrix cannot be consistent giving no 
solution or infinite solutions (open-loop solutions) 
cases. Thus, singularity of matrices occurring both 
naturally and those arising during computations by 
truncation of numbers to a finite number of digits, 
are a hurdle towards a_ successful analytical 
solution. 

A singular or degenerate matrix can appear to be 
non-singular on cursory examination. The 
prechecking by computation of determinant|A| + 


0 can immediately resolve the matter by as we have 
seen (section on Cramer’s rule) this is 
computationally an inefficient method of solving 
matrices as observed by Proskuryakov, I.V. (2011). 
A matrix is [A]is non-singular on the condition that 
the diagonal elements are nonzero stated by Golub, 
G. H., Van Loan, C. F. (1996). A condition which can 
be easily checked by diagonal element of a;j; € 
[A] i (row) = j (column) and a;,; # 0 condition in 
loop. This cursory evaluation will later fail as 
singularity may arise due to loss of precision 
especially while rounding errors due to truncation 
of digits or bits in floating point operation, such 
matrices are said to be ill-conditioned. Prechecking 
by finding the Determinant of a large array is 
wasteful and not efficient for large matrices as 
described by Xian-Da Zhang (2017). Since 
definition of limit of singularity in calculating 
matrices is defined ne a — 0, we may immediately 


note that the magnitude of the determinant is the 
simple quantification of an ill-conditioned a matrix. 
lim|A| > 0 
n 


A more precise, Perturbation theory-based formula 
for condition of a matrix (ill-behaved, singularity) is 
the computation of formula described by Sewell, G. 
(2017). 


conditionof [A] = |lAl|||A~*|| > 1 
This loops us back to the paradox of Matrix LU 
decomposition, which is illustrated in flowchart 
Fig.3. 
A simple example of ill-conditioned matrix which 
produce singularity in computations, benign 
looking Hilbertian matrix as demarcated by Costa, 
R. (2020). Defined by this equation 
Hilbertianmatrices[H] = aay forts = 
1,2,3,4 soonN forlargeN, 
12 3 
2 3 i and lim|H| — 0 (Sewell, G. 2017). 
3.4 °5 


a smaller example is 


Lotkin’s matrices which are identical to the 
Hilbertian matrix, first row is completed with 1 


only (Sewell, G. 2017). 
114 
[L] = E 3 i 
3.4 5 
Using floating-point arithmetic of computer system 
when an operation of summation of n numbers is 
performed, each number contributing to an 
infinitesimally small number 6d, -:-dd, computed 
the exact sum of the same n numbers perturbed by 
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the computed analytical closed-loop solution differs 
from the real solution obtained by hand. This 
happens because 
lim107-" = 0 «. 


n-0 
degenerate and solutions have become open loop in 


midst of computation. 


infinitesimally relative errors. The errors add up in 
computation like this. 


> X,(Computedwitherrors) 


= X,(1 + 6d,)+....+X,(1 + 6d,) 
Loss of Precision 
Since floating point operations on computer 
processor are only performed to a limited precision 


Matrix has become Singular / 


Sucocssful LU Factori- 
zaboe needs chat Matrix 
[A| de man-singulur 


To find _|A) is not vero in an efficient 
manner needs Successful LU Factori- 
zation of [|A| Direct, classroom- 
hand solution not efficient for 
large matrix A as described) 


Determinant |A{ ma- 
trixes [L| [U] matri- 
ces (= diagonal 
product) 


Interchange rows moving highest 
values mn cuch row to diauonal so 
that sinwularities while competing 
are avoeded by partial prvoting), 
use Permutation mutrix [P} 


Is the [A] Matrix singu> 
lar on pecliminury in 


spection’? 
NO 


Check for condition 
Oy, has some arbitrary 


i 


Proceed with LU Fae- 
sorizaton of the Matrix 


singularity Threshold value 


I) Matrix [A] has bocome singu__ 
tar due to loss of Precision 


No solution [A], equa- 
tion is unsolvable 
NOW duc to Loss of 
Precision 


Fig4. The algorithm of avoiding Singular matrices in large matrix computations and the paradox of 
LU factorization and solving the Determinant of a large matrix. This paradox is narrated thus, the 
matrix can be factorized, efficiently if determinant|A| + 0 but to find the determinant|A| we need 
to factorize the matrix successfully and find L and U matrices (Xian-Da Zhang 2017). The 
singularity threshold of Apache Java Math library (commons-math 3-3.6) is arbitrarily set to 
singularity threshold 10-1! ~ 0 is effectively zero, this I tested on Mobile CPU, ART Android 


Runtime (above KitKat) and found that actual truncation errors began at orders > 107 


308 only but 


this might be a feature of the ART machine of the benchmarked tested mobile. 


LU factorization by partial pivoting using a 
Permutation Matrix [P] 


Pivot is the first and the left most nonzero element 
in a row of a matrix, beginning at which is the 
starting point of elimination by Gauss method. We 
find that the first pivot element a1,1 =0 or, 
lim a1, ~ 0, we need to Pivot or rearrange the 


rows so that the row a, # 0, this process is called 
Partial Pivoting sated by Kaw, A.K. (2008). 

Interchange of both rows and columns is complete 
pivoting. Complete pivoting is more reliable in 
avoidance of matrix singularity, involves testing 
each element in row and column. Conversely, this 
places considerable computational overhead, 
destroys the structure of a matrix and cannot be 


employed on blocked matrices and so is not 
considered an ideal method stated by Sewell, G. 
(2017). 

Partial Pivoting for size is agreed best tradeoff 
method in terms of accuracy, algorithm stability, 
computational efficiency works well real-life 
situations. The concept of pivoting is simple, but it 
is not computationally trivial matter to decide 
which element to use as pivot. A row is searched, 
and number of largest magnitudes, i-e., modulus in 
case negative number is found and moved to the 
location (1,1) in the matrix. This step is repeated 
successfully in each row as described by Malyshev, 
A. N. (1995) and Golub, G.H.; Meurant G.A. (2006). 
Other pivoting strategies are Scaled pivoting by 
calculating ratios and Rook Pivoting. 
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The existence of such a row is attested by that fact 
that determinant of [A], |A| # 0. This is possible by 
left product of A by a Permutation matrix [P], 
whose properties we already summarized. The 
proposed analytical, solution for the Linear 
equation of N size is as follows. 

[P] x [A] = [L] x [U] 

Permutation matrix is a simple binary matrix (of 
only a=1 and Os) to allow switching of rows, and 
please note that inverse matrix of [P]~tis same as 
[P]. 

[P] = inverseof [P| 

Or 

[A] = [L] x [U] 

Analytical closed-loop Solutions of system of Linear 
equations using LU factorization is by following 
steps 

[A] x [X] = [C] 

[L] x [U] x [X] = [P]~* x [C] = [P] x [C] 
Substitute a column matrix [Y] = [U] x [X] 


Then forward solve, called forward sweep, 
substituting [L] x [Y] =[P] x [C] for values in 
matrix [Y] 


Then backward solve, called backward sweep 
substituting matrix [Y] in [U] x [X]=[Y] or 
finally 
[x] = [U]* x [Y] 
The forward sweep places an overhead of N(N — 
1) floating-point operations and backward sweep 
places N(NW — 1)+WN floating point operations, 
compared this to the actual LU decomposition step 
which is O(N?) as clarified by Meurant G.A. 
(2006) 
Choice of the available library Apache Math 
Common 
Apache Math Common (commons-math3 version 
3.6) is a Java language library of lightweight, self- 
contained, common mathematical and _ statistics 
libraries for solving common computer science 
problems, available from the website 
https://commons.apache.org/proper/commons- 
math/. 
Most of the algorithms used in this package are 
optimized and best practice to utilize least amount 
of system resources in Computation and provide 
reliable results on desktop and server versions of 
Java. We could not find any previous reference or 
experience to porting this desktop Java library to 
the ART Android Run Time Processor (C.P.U). 


Issues in actual implementation of the LU 
Decomposition 


I faced the following issues / problems in actual 
implementation of the algorithm 

Limited memory RAM and CPU computation power 
of a cell phone to perform Matrix computations 
(Supercomputers were used) 

Create an error free GUI interface (also use a lot of 
9-patches as background for Input dialogs) 
Validation of user input 

Check is the Matrix is square by checking number 
of row equal to each column of the Matrix (cannot 
solve a ragged 2D array) 

Check is the system of linear equations has an 
analytical, closed-loop solution (reached within 
finite number of steps) 

Check is the Matrix is singular on cursory 
evaluation 

Feasibility of porting the Apache Java library to 
Android (described here first time there is no 
documentation extant) 

Choosing Partial pivoting row to avoid singularity 
of matrix in midst of calculations due to retention 
only a finite number of digits (or bits) in floating 
point calculation using singularity threshold 
1071! = O is effectively zero. 

Floating point operation without leaking memory 


Accurate output 


Complete Android Test and Lint error free, push 
App to debug stage 

Use Object Oriented programming (no go to 
jumps!) and reuse cases 

Deal with issue of obfuscation of Java classes by 
ProGuard (dealt with successfully for first time in 
this paper), push Java compilation to App release 
stage 

Successful error free compilation on Android studio 
Document LU factorization in terms of speed and 
accuracy 

Modelling Matrix singularities due to truncation in 
Android Run Time using Poisson Distribution 
Poisson distribution is used in analysis of events 
which are sporadic and have a rare likelihood of 
occurrence. A simple rule of thumb for choosing 
Poisson distribution as a decent approximation to 
Binomial distribution if N (number of variables or 
computations both criteria use in this study) is 
larger than > 20 and NXP <«10o0reven<1 
which is true in our study. Poisson distribution is a 
limiting case of Gaussian distribution when the 
probability of a rare event P, LimP-0O, is 
diminishingly small and the LimQ - 1, Q is large, 
limiting to 1, where is P is probability of rare event 
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happening and Q probability of event not occurring. 
A specimen of rare event P is Truncation errors 
leading to singularities despite using partial 
pivoting procedure. Poisson probability of 
obtaining N positive outcomes is given by 

i xe? 
P(NraregventErrors) = NI 
Where py is population mean or also same as 
variance related this way o7=NxPx 
QandLimQ —1,LimP-—0 then we can deduce 
thato7=N xP = p. 


Materials 


Graphical User Interface development (G.U.I) on 
Android Studio development process 

The program for implementation of LU 
Decomposition with partial pivoting was done on 
Android Java version 8, Android Studio version 
Gradle Version 1.5.0 and JetBrains Editor. Release 
compilation was done targeting Android Software 
development kit Version 23.0, though application 
runs fine on any Android mobile phone version 
above 9 (Google Developers, 2019). This was 
achieved by including appcompat library version 
v7-23.1.1, to ensure backward compatibility with 
older versions of Android phones (Android 
Developers, 2021). 

ProGuard is a Java optimizer for Android apps 
which makes execution 20% faster and compresses 
app up to 90%. ProGuard also protects against 
decompilation and reverse engineering of the 
application. ProGuard works for both Java and 
Kotlin. ProGuard's Gradle plugin was applied via 
the command in the file build.gradle, using 
command, apply ‘plugin: 
com.guardsquare.proguard'. ProGuard version was 
applied to obfuscate the classes which caused the 
loss of classes leading to crashing the app on beta 
testing. To avoid obfuscation of essential classes in 
Math 3, Linear Algebra library for successful 
operation, the following lines were added to the 
text file which contains customized rules for 
application of ProGuard, proguard-rules.pro 
(www.guardsquare.com, n.d.). 

-keep class org. apache. commons. math3.linear. ** 
-keep interface org. apache. commons. 
math3.linear. ** 

-dontwarn org. apache. commons. math3.geometry. 


**K 


The release version of application was also Zip 
Aligned by the command, zipAlignEnabled true and 
Gradle command ‘gradlewassembleRelease’ was 
used to compile the release candidate application. 
(www.guardsquare.com, n.d.) 

A 9-Patch image is a background image which can 
dynamically accommodate with contents of a view 
like text or image without pixelating or blurring. 9- 
patches are very useful in highlighting the Edit 
text/ Input fields for user input. 9-pathches were 
used extensively for building the GUI of the 
application. These are capable to expanding to 
accommodate increasing amounts of numerical 
input. left and top are for scaling, and the right and 
bottom are for content. The left and top margins of 
a rectangular image will stretch to fit, and the 
contents will be shown padded by the black lines 
along the right and bottom margins. Android studio 
provides a Java Standalone desktop application for 
creating and editing 9-patches. 

Android Application was successfully complied, and 
release version produced and uploaded to the 
GitHub.com site available at also at online research 
repository Figshare.com as for general repository, 
links are available in “Data availability statement” 
section at end of this paper. 


Operation of the Graphical User Interface 
development (G.U.I) 


The G.U.I is an Android drawer layout (Android 
Library recyclerview-v7-23.1.1) with a floating 
action button titled “RUN”. The matrices are input 
by the user using the soft keyboard that pops up on 
touching the green area which is the Matrix [A]. The 
Matrix [C]or the Column vectors on right side of 
matrix equation [A][X] = [C] are fed on the blue 
area. Matrix columns elements are separated by 
single Space key press and rows elements are 
separated by an Enter (J) key press. The linear 
equation of N size can then be solved by pressing 
the “RUN” button in right lower hand corner. The 
program then proceeds to convert the String array 
to into a 2-dimensional array or matrix. After 
cursory validation for Singulate matrix, the 
analytical solution of is depicted in the ListView 
along with Determinant of a N size matrix equation. 
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Graphical User Interface G.UI 


Launcher icon on 
Home screen 


Execute the 
matrix computation 
by “RUN” 
Floating action 
button 


Listview that 
shows unique, non-trivial 
solutions for any matrix [A] 


Fig5. The Graphical user interface G.U.I of the release app that computes the unique closed-loop, 
analytical solutions to dense linear, algebraic equations represented by matrix equation [A][X] = 
[C], with N unknown variables by LU Decomposition with partial pivoting algorithm on an Android 
mobile. The blue matrix and green matrix are 9-patches which will dynamically resize to enclose 
the user input into the matrices. The green area is the left-hand side of the equation matrix [A] and 
the blue area is the right hand side matrix [C].Matrix [X] of the unknowns (41, %2,%3,""" Xv = 
N variables)are automatically inferred from the N= size of the rows or columns of square matrix 
provided by the user. Additional padding with zeros is needed if the coefficient of any linear 
equation is zero like so 0 x x, so that the matrix [A]is converted into a square with rows = 
columns. This app is hosted on the GitHub platform for free download and use 
(https: //github.com /Kephalian/Android-Dense-Linear-Equation-Solver). 


Observations 


Benchmarking and Beta testing of the Application 
Since matrix computation are most resource 
consuming these were done on Supercomputers 
and Main frames few decades ago, now these same 
number-crunching, intensive matrix is capable of 
being executed on an ordinary, middle range 
Samsung A30 mobile phone [26]. 

The objective of this benchmarking was _ to 
demonstrate the robustness, validity, accuracy, and 
efficiency of matrix processing on the Android CPU. 
1.Test the robustness of LU Decomposition 


algorithms with partial pivoting on a mobile 
processor 

2. To test accuracy by computing solutions with 
extremely small values in matrix [A] and 
occurrence computational singularity arising out of 
truncation digits and to model the outcomes on a 
Poisson distribution mass density function 

3.To analyze the Execution speed of the program 
and correlating with Number of equations (N), 
Number of floating-point operations and Size of the 
matrix [A] (N x N). 
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Fig6. Modelling the incidence of Singularities in midst of matrix computation by feeding 
progressively smaller values in the Matrix LU decomposition algorithm. The density of successful 
analytical closed-loop solution (Poisson variable Q—1) is plotted, and outlier Poisson event P is 
shown. Since the occurrence of singularity is a stochastic, discrete, time event it is modelled best 
by Poisson distribution. Thermal 3D surface plot of the density shows increasing events with time, 
depicted by upslope, and increased red colour (like increasing temperature) confirming that these 

errors are a type of time-dependent Poisson events. 


Table1. Frequencies of Matrix singularity observed in 2 Poisson event: Poisson event 1, Q Limit >1 
is successful matrix computation and Poisson event 2, Limit P— 0 is failure to compute due to 
accumulation of rounding errors caused machine truncation of digits in floating-point operation. 
Error rate in the machine is calculated to be 51 out of 250 x 1000 variables = 2.04 x 10-°per N 
(variable) and gross error rate per computation is calculated by 1000 x N? = 1000 x 250? = 
15,625,000, 000 fops hence 3.264 x 10~°per floating point operations. 


Test Benchmark test of computation | Incidence | Computed P value following 
Numbe of error | Poisson event P of error, N? 
r or fops performed 1000 times 
Singularit 
y out of 
1000 
runs for 
unique 
analytical 
solutions 
N= 10 Equation at limit 107308 4 4x 107° 
N=10 Equation at limit10~27° 3 3x 10°° 
N=10 Equation at limit 1077°° 6 6x 10° 
N=10 Equation at limit 107148 3 a6 10" 
N=30 Equation at limit 10~1°° 9 3.33 x 1077 
N=30 Equation at limit 10~°° 3 LA 10-7 
N=30 Equation at limit 10778 5 1.852 x 10-7 
N=30 Equation at limit 10~7° 3 11d 10>" 
N=30 Equation at limit 10~7° 4 1.48 x 1077 
Singularity Benchmark N=30 | 3 cs te a 
Equation at limit 10~® 
N=30 Equation at limit 107? 8 2.96 x 1077 
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Total 11 tests X 1000 runs each | )}.X,(1 + | Porror error) 
X total N=250 X 1000 6d,)+....4 = 2.4 
6d,)=51 | x 10-*pervariablecomputed, 
Porror (2 errors) is 
= 2.08 
Model x 10-8 perunknownvariablex| 
ling Porror(1 error) = 3.26 X 
Singul 10-°Per floating-point 
ar operation 
matri Porror (2 errors) is = 5.34 x 
x 10718 per fops 
incide 
nce in computation of non-trivial solutions by 
inputting extremely small, randomized values in were divided serially by 


elements matrix [A] 

Random number generator was used to populate 
the matrix of [A], and[C]( [X] is inferred from size 
of N). The number of variables in the test were 
initially linear algebraic equations of 30 variables 
(square matrix size [A]isN?=300), later reduced 
equations of 10 variables (Matrix size 100). Then 
these values 


1S 4 -— Loganthmic Posson distribution plot P{Simguiarities 


— Loganthic Posse 


10 A010" a A ee 


Or analy" 

, in an attempt to stimulate computational 
singularity by truncation of errors. The results are 
shown in table below. Beyond 1073°% we 
encountered ART Android Run Time machine error 
due to underflow and analytical solutions are not 
possible on the tested model of mobile phone. 


v.s {1 Adtnals tor N Venables 


mn destnbution plot PiSemquiantes) vs (1k) traits FOPs Operations 


Fig7. Mathematical Modelling of the ART CPU truncation errors leading to Matrix singularity based 
on Poisson distribution of K trails showing asymmetrical curve, skewed towards the right, or 
positively skewed as probability density is nearly zero (as P « Q,u = 0% = N x P), Logarithmic 
probabilities P of K trails were plotted against 1/k as P is very small. This skewed curve has 


Kurtosis of — = 70.014 and 17503.50. 


VE 


Modelling based on the data available P.,,9-0TL = 
2.04 x 10°°%per N (variable) and P,,,9,0TL = 
3.264 x per 10~°fop in Poisson model we can see 
that. 


po sce? 
P(Nyarekventkrrors) = ae 
so Porror(1 error) = 2.04 x 


10-*pervariablecomputed, Poyyoy(2 errors) is = 


2.081 x 10-8per unknown x solved, in other terms 
of Floating-point operations per _ operation 
performed on ART machine P,,,;9,(1 error) = 
3.26 x 10-°Per floating point operation 
performed on Android Run _ time 
machine,Poryor(2 errors) is = 5.327x 10718 per 
fops. Indeed, this shows’ the advanced 
computational abilities of modest, cell phone 
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processor and the robustness of the algorithm 
employed. 

Analyses of the Execution speed of the program and 
correlating with Number of equations (NV), Number 
of floating-point operations and Size of the matrix 
[A] (N x N). 

The execution speeds were tested for benchmark 


by increasing the number of variables solved from 
N = 2 to 60 variables, computation N3 = 8 fops to 
2,16,000 fops in a linear algebraic equation run X 
1000 times (8000 fops to maximum mind-boggling 
2,16,000,000 fops). The following are descriptive 
data regarding the benchmark tests. 


Table2. Average execution times in psecondsfor solving dense, linear equations of [A] on Android 
mobile phone Samsung A30. The number of runs was 100 and the number of variables varied from 
N= 2 unknowns, needing 8000 fops to N=60 unknowns needing 2,16,000,000 fops or accurate 


steps. 

Serial System of 2 Linear | System of 4 Linear | System of 10 Linear 
No. Equations LU | Equations LU | Equations LU 
Decomposition Decomposition Decomposition 
microseconds 107° microseconds 107° microseconds 107° 

Mean 0.0698 0.173 1.004 
Executi pseconds pseconds puseconds 

on 

Time 

Std. 0.610 0.694 1.51 

error 

mean 

Maximu 50.6 127 591 

m pseconds puseconds piseconds 

Minimu 246 375 1315 

m pseconds puseconds puseconds 
System of 40 Linear | System of 50 Linear | System of 60 Linear 
Equations LU | Equations LU | Equations LU 
Decomposition Decomposition Decomposition 
microseconds 107° microseconds 107° microseconds 107° 

Mean 29.213 54.164 88.524 

Executi pseconds pseconds pseconds 

on 

Time 

Std. 173 431 604 

error 

mean 

Maximu 8105 12567 17969 

m pseconds puseconds puseconds 

Minimu 43968 227375 99430 

m pseconds pseconds pseconds 


The Karl Pearson's r correlation coefficient (r) 
between the variables Mean execution times and 
Number of variables N was calculated to be + 0.938 
(*** p < .001), in addition Number of Floating-Point 
Operations Performed in calculation N?and 
Number of Linear Equations with N Unknowns was 
calculated to be + 0.924 (*** p < .001), 

As logically as expected to be a positive, linear 


correlation, which means that the execution time 
increases linearly with the increasing complexity of 
calculations required to solve the matrix equations. 


Conclusion 


This paper is a distillation of my experience in 
building and testing an Android application which 
provided closed-loop analytical solution to dense 
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linear equations which contain any arbitrary 
number of variables for the purpose of modelling 
and network analysis in healthcare and 
neuroscience research. The paper begins by 
introducing the concept of network analysis in 
modelling for healthcare and neuroscience 
research. A simple example of network modelling is 
the logistics of delivering vaccine like COVID19 
vaccine from a company to Primary health center 
while maintaining cold storage of the vaccine. Next, 
a brief algebraic background on Matrix equation 
methods of solutions, the paper describes the 
paradigm of LU decomposition with partial pivoting 
including its triumphs and drawbacks. The 
implementation of LU factorization of matrices on 
computational machines with truncation of digits is 
a hurdle in successful and resource efficient 
factorization. The paper describes the process in 
development of the Android application and 
concerns in building and compiling the release 
application to Android machine code. Testing the 
robustness, efficiency, and reliability of this 
implementation on an Android Samsung mid-range 
mobile phone is described next. Computational 
singularity event or error of Matrix computations 
encountered in solving matrix algebraic equations 
are rare, independent, discrete, stochastic events 
and are mathematically modelled on Poisson 
Probability mass function. The calculated Poisson 
modelling of probability of Stochastic, Singularity 
event is = 3.26 x 10~°per floating point operation 
a number derived after running 15,625,000,000 
floating-point operation runs. The mean execution 
time was 88.534 seconds , for solving matrix 
equation [A]withN=60, 60 variables and 3600 rows 
and columns, performing N?=216000 
computations. As anticipated, the CPU execution 
times show statistically significant large, positive 
correlation (Karl Pearson's r correlation coefficient 
(r) is + 0.938, p value < .001) with increasing 
complexity of equations (N) and sizes of matrices 
(N*). 


Data Availability Statement 


Android Application was successfully complied, and 
release version produced and uploaded to the 
GitHub site for free download, and available at 
online research repository Figshare.com 
(https://figshare.com/) as per Springer 
publications guidelines for general repository at the 
following site: 
https://github.com/kephalian/Android-Dense- 


Linear-Equation-Solver 

The compiled Android application, datasets, 
statistical analysis files, extra graphs, Java, and 
Python sources generated during the current study 
are available in the _ (https://figshare.com/) 
repository: 
https://figshare.com/projects/Android_Applicatio 
n_Dense_Linear_Equation_Solver_using_LU_decomp 
osition_and_partial_pivoting_for_mobile_phones/14 
3931 
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